---
title: "P4D Initial Analysis"
author: "Jonathan Pinckney"
date: "1/13/2021"
output: html_document
---

## Introduction

The following RMarkdown will replicate all the numbers and figures included in "Not Just Who but How: Getting from Nonviolent Action to Democracy Through Inclusive Dialogue and Negotiation Processes." To ensure the code runs correctly, all the packages in the set-up chunk should be installed, and the "data_prep.R" script should be in the same working directory as this markdown file. 

You may run the individual chunks of the markdown in order to reproduce the output, or knit the entire markdown to produce an html document with all the results.



```{r setup, include=FALSE}

library(tidyverse)
library(rnaturalearth) # Package for creating maps
library(sf) # Package for spatial analysis
library(lubridate) # For working with dates
library(margins)
library(broom)

# Additional Packages Needed for Data Prep Script

# install.packages("readxl")
# install.packages("countrycode")
# install.packages("devtools")
# devtools::install_github("vdeminstitute/vdemdata")

source("data_prep.R") # Run data prep script

# Get Color Hex Codes for USIP color scheme

usip.yellow <- "#C19F53"
usip.blue <-  "#276383"
usip.dark.blue <- "#183B4B"
usip.gray <- "#949599"


usip.red.brown <- "#C18467"
usip.other.red.brown <- "#856466"
usip.gray.green <- "#68836E"
usip.light.gray <- "#ECE7DA"
usip.dark.teal <- "#36868D"
usip.light.teal <- "#98C3C8"

```

## Summary Statistics

### Map of Civil Resistance Transitions

```{r}
crts.for.map <- crts %>% 
  mutate(cowcode = case_when(cowcode == 678 ~ 679, # Change a few COW codes to match ISO codes 
                             cowcode == 265 ~ 255,
                             cowcode == 347 ~ 345,
                             TRUE ~ cowcode),
         iso = countrycode::countrycode(cowcode,'cown','iso3c',warn = T), # Add ISO codes
         crt = 1) %>% # Create indicator of a CRT
  distinct(iso,.keep_all = T) %>% # Remove duplicate lines of countries that have had more than one CRT
  mutate(iso = if_else(cowcode == 345,"SRB",iso))


mapdata <- ne_countries(returnclass = "sf") %>%  # Pull in simple features data frame of world map
  left_join(.,crts.for.map,by = c("iso_a3" = "iso")) %>% # Join to CRT list
  mutate(crt = if_else(is.na(crt),0,crt)) %>% # Make indicator of CRT == 0 if no CRT in that country
  filter(admin != "Antarctica") # Remove Antarctica from map


crtmap.usip <- ggplot() + geom_sf(data = mapdata,aes(fill = as.factor(crt)),color = "gray37") +
  scale_fill_manual(name = NULL,
                  breaks = c(0,1),
                  values = c(usip.gray,usip.yellow),
                  labels = c("No NVA Transition","NVA Transition")
                  ) +
  coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") +
  theme_bw() +
  theme(panel.background = element_rect(fill = "lightblue"),
        panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        text = element_text(size = 36),
        legend.position = "bottom")

crtmap.usip

mean(p4d_data$index)

# ggsave("figures/usip_map.jpg",plot = crtmap.usip,height = 10,width = 20,units = "in",dpi = 800)
```


## Duration and Scope

### Duration

```{r}

mean(p4d_data$duration) # Mean is 250 days

median(p4d_data$duration) # Mean is 250 days

mean(p4d_data$agreement_bin) # Roughly 71%

ggplot(filter(p4d_data,duration < 1000),aes(duration)) + 
  geom_histogram(color = "black", bins = 50, fill = usip.blue) +
  labs(y = "Number of DNPs",
       x = "Length in Days",
       caption = "6 DNPs over 1,000 days excluded for ease of visualization") +
  theme_bw()

```


### Distribution of Scope

```{r}

sum(p4d_data$scope_new %in% c("revolutionary","constitutional"))/nrow(p4d_data)

ggplot(filter(p4d_data, !(is.na(scope_new))),aes(factor(scope_new, levels = c("revolutionary","constitutional","policy","no mandate")))) + 
  geom_bar(fill = usip.blue) +
  scale_x_discrete(labels = c("Revolutionary","Constitutional","Policy","No Mandate \n For Change")) +
  theme_bw() +
  labs(y = "Number of DNPs", x = "Scope")


```

## Who participates in DNPs?

```{r}
# Participation of organized groups

p4d_data %>% 
  summarize_at(vars(groups.oppo:groups.gov), ~ sum(.)) %>% 
  pivot_longer(cols = groups.oppo:groups.gov) %>% 
  mutate(value = value/nrow(p4d_data),
         name = c("Opposition \n Parties","CSOs","Religious \n Groups","Ethnic \n Organizations","Military","Pro-Gov \n Parties","Nonstate \n Armed Groups","Government")) %>% 
  ggplot(aes(reorder(name,desc(value)),value)) + geom_col(fill = usip.blue) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "Participating Groups",y = "Percent of DNPs Included") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45,hjust = 1))

# Participation of Women

mean(p4d_data$women_bin,na.rm = T)

table(p4d_data$women)

```

## Selection Mechanisms

```{r}
p4d_data %>% 
  summarize_at(vars(selection.govt:selection.appointed), ~ sum(.,na.rm= T)) %>% 
  pivot_longer(cols = selection.govt:selection.appointed) %>% 
  mutate(value = value/nrow(p4d_data))
```

## Balance of Power, decision-making mechanism, and agreement characteristics

```{r}
sum(p4d_data$domination > 2,na.rm = T)/nrow(p4d_data) # 70% have at least balance between old and new

sum(p4d_data$domination == 5,na.rm = T)

table(p4d_data$decisions)

mean(p4d_data$provisions.key.issues,na.rm = T)


```

### Prevalence of Mediation

```{r}
mean(p4d_data$int.mediation.bin,na.rm = T)

ggplot(p4d_data,aes(int.mediation.bin)) + geom_histogram(stat = "count")

full.dataset %>%  # Prevalence of International Mediation Alone
  group_by(decade,any.intl.mediation,any.dom.mediation) %>% 
  summarize(count = n()) %>% 
  ggplot(aes(x = decade,y = count,fill = factor(any.intl.mediation,levels = c("0","1")))) + geom_col() +
  scale_fill_discrete(name = NULL,labels = c("No Mediation","Mediation")) +
  labs(x = "Decade",y = "Number of Transitions") +
  theme_bw()


full.dataset %>%  # Prevalence of Intl and Domestic Mediation
  group_by(decade,any.intl.mediation,any.dom.mediation) %>% 
  summarize(count = n()) %>% 
  mutate(Mediation = case_when(any.intl.mediation == 1 & any.dom.mediation == 1 ~ "Both",
                                   any.intl.mediation == 1 & any.dom.mediation == 0 ~ "Intl. Only",
                                   any.intl.mediation == 0 & any.dom.mediation == 1 ~ "Dom. Only",
                                   TRUE ~ "None")) %>% 
  ggplot(aes(x = decade,y = count,fill = Mediation)) + geom_col() +
  # scale_fill_discrete(name = NULL,labels = c("No Mediation","Mediation")) +
  labs(x = "Decade",y = "Number of Transitions") +
  scale_fill_manual(breaks = c("Dom. Only","Intl. Only","Both","None"),values = c(usip.light.teal,usip.blue,usip.yellow,usip.gray)) +
  theme_bw() 

mediation.plot <- full.dataset %>%  # Prevalence of Intl and Domestic Mediation
  group_by(decade,any.intl.mediation,any.dom.mediation) %>% 
  summarize(count = n()) %>% 
  mutate(Mediation = case_when(any.intl.mediation == 1 & any.dom.mediation == 1 ~ "Both",
                                   any.intl.mediation == 1 & any.dom.mediation == 0 ~ "Intl. Only",
                                   any.intl.mediation == 0 & any.dom.mediation == 1 ~ "Dom. Only",
                                   TRUE ~ "None")) %>% 
  ggplot(aes(x = decade,y = count,fill = Mediation)) + geom_col() +
  # scale_fill_discrete(name = NULL,labels = c("No Mediation","Mediation")) +
  labs(x = "Decade",y = "Number of Transitions") +
  scale_fill_manual(breaks = c("Dom. Only","Intl. Only","Both","None"),values = c(usip.light.teal,usip.blue,usip.yellow,usip.gray)) +
  theme_bw() +
  theme(legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(fill = "transparent",color = "transparent"),
        plot.background = element_rect(fill = "transparent",color = "transparent"),
        text = element_text(color = "white",size = 18),
        axis.text = element_text(color = "white",size = 14))

ggsave("mediation_plot.png",mediation.plot,bg = "transparent",width = 9,height = 6,units = "in")



```

## Democratization

### Inclusion Index Models

```{r}
index.mod1 <- lm(v2x_polyarchy_lead5 ~ index + v2x_polyarchy_lag1 + e_migdppcln + region_polyarchy,data = full.dataset)
index.mod2 <- lm(v2x_delibdem_lead5 ~ index + v2x_polyarchy_lag1 + e_migdppcln + region_polyarchy,data = full.dataset)

summary(index.mod1)

fig.data <- tibble(index = 2:6,
                   v2x_polyarchy_lag1 = mean(full.dataset$v2x_polyarchy_lag1, na.rm = T),
                   e_migdppcln = mean(full.dataset$e_migdppcln,na.rm = T),
                   region_polyarchy = mean(full.dataset$region_polyarchy,na.rm = T)) %>% 
  mutate(poly.predicted = predict(index.mod1,newdata = .),
         poly.predicted.lwr = predict(index.mod1,newdata = .,interval = "confidence")[,2],
         poly.predicted.upr = predict(index.mod1,newdata = .,interval = "confidence")[,3],
         delibdem.predicted = predict(index.mod2,newdata = .),
         delibdem.predicted.lwr = predict(index.mod2,newdata = .,interval = "confidence")[,2],
         delibdem.predicted.upr = predict(index.mod2,newdata = .,interval = "confidence")[,3]
         )

ggplot(fig.data,aes(x = index,y = poly.predicted,ymin = poly.predicted.lwr,ymax = poly.predicted.upr)) +
  geom_ribbon(fill = "lightgray") +
  geom_line() +
  scale_y_continuous(name = "Predicted Polyarchy",
                     limits = c(0,0.7)) +
  scale_x_continuous(name = "Level of Inclusion",
                     breaks = c(2,4,6),
                     labels = c("Low","Medium","High")) +
  theme_bw()

ggplot(fig.data,aes(x = index,y = delibdem.predicted,ymin = delibdem.predicted.lwr,ymax = delibdem.predicted.upr)) +
  geom_ribbon(fill = "lightgray") +
  geom_line() +
  scale_y_continuous(name = "Predicted Democracy",
                     limits = c(0,0.6)) +
  scale_x_continuous(name = "Level of Inclusion",
                     breaks = c(2,4,6),
                     labels = c("Low","Medium","High")) +
  theme_bw() +
  theme(axis.text = element_text(size = 14),
        text = element_text(size = 18))


# Fig Data Bar Chart (Index and DNP only)
dnp.fig.data.bar.chart <- dnp.fig.data %>% 
  ggplot(aes(dnp.binary,delibdem.predicted,fill = as.factor(dnp.binary))) + 
  geom_col(width = 0.7) +
  scale_fill_manual(values = c(usip.blue,usip.yellow))+
  scale_x_continuous(limits = c(-0.5,1.5),breaks = c(0,1),labels = c("No DNPs","DNPs")) +
  labs(y = "Deliberative Democracy Score",x = "") +
  theme_bw() +
  theme(legend.position = "none",
        text = element_text(size = 18))

dnp.fig.data.bar.chart

index.fig.data.bar.chart <- index.fig.data %>% 
  filter(index %in% c(2,4,6)) %>% 
  ggplot(aes(index,delibdem.predicted,fill = as.factor(index))) + 
  geom_col() +
  scale_fill_manual(values = c(usip.dark.blue,usip.blue,usip.yellow))+
  scale_x_continuous(breaks = c(2,4,6),labels = c("Low","Medium","High")) +
  labs(y = "Deliberative Democracy Score",x = "Level of Inclusion") +
  theme_bw() +
  theme(legend.position = "none",
        text = element_text(size = 18))
index.fig.data.bar.chart

```



### Disaggregated Models

```{r}
mod1 <- lm(v2x_polyarchy_lead5 ~ dnp.binary + domination.mean + any.women.partip + any.cso.partip + any.agreement +
           v2x_polyarchy_lag1 + e_migdppcln + region_polyarchy, data = full.dataset)

summary(mod1)

coef.order <- c("e_migdppcln","region_polyarchy","v2x_polyarchy_lag1","domination.mean","I(domination.mean^2)","any.cso.partip","any.women.partip","any.agreement","dnp.binary")

# Coefficient Plot

tidy(mod1,conf.int = T) %>%
  filter(term != "(Intercept)") %>% 
  mutate(significant = if_else(p.value < 0.05,"Sig","Nonsig"),
         term = factor(term,levels = coef.order)) %>% 
  ggplot(aes(term,estimate,color = significant)) + geom_point() +
  geom_errorbar(aes(ymin = conf.low,ymax = conf.high,color = significant),width = 0.2) +
  geom_hline(yintercept = 0,linetype = "dashed") +
  coord_flip() +
  scale_color_manual(values = c("black","red"),labels = c("No","Yes")) +
  theme_bw() +
  labs(x = "",y = "")
    
```
## Women's Participation Predicted Value 

```{r}

tibble(any.women.partip = 0:1,
                   dnp.binary = 1,
                   e_migdppcln = mean(full.dataset$e_migdppcln,na.rm = T),
                   region_polyarchy = mean(full.dataset$region_polyarchy,na.rm = T),
                   v2x_polyarchy_lag1 = mean(full.dataset$v2x_polyarchy_lag1,na.rm = T),
                   domination.mean = mean(full.dataset$domination.mean,na.rm = T),
                   any.cso.partip = 1,
                   any.agreement = 1) %>% 
  mutate(poly.predicted = predict(mod1,newdata = .)) %>% 
  ggplot(aes(any.women.partip,poly.predicted,fill = as.factor(any.women.partip))) + 
  geom_col() +
  scale_fill_manual(values = c(usip.blue,usip.yellow))+
  scale_x_continuous(breaks = c(0,1),labels = c("No Women's \nParticipation","Women's \nParticipation")) +
  labs(y = "Polyarchy Score",x = "") +
  theme_bw() +
  theme(legend.position = "none")


```
## Case Study Countries

```{r}

case.studies <- full.dataset %>% 
  filter(Startyear %in% c(2014,2011) & Country %in% c("Tunisia","Egypt","Ukraine")) %>% 
  select(Country,v2x_polyarchy_lead5,v2x_delibdem_lead5) %>% 
  pivot_longer(-Country) %>% 
  ggplot(aes(x = factor(Country, levels = c("Tunisia","Ukraine","Egypt")) ,y = value,fill = name)) +
  geom_col(position = position_dodge()) +
  labs(y = "Democracy Score",x = "") +
  scale_fill_manual(values = c(usip.blue,usip.yellow),labels = c("Deliberative\nDemocracy","Polyarchy"),name = "") +
  theme_bw() +
  theme()

case.studies
```

